COMMENT ⊗   VALID 00006 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002		SUBROUTINE LPC(AIFFY,SPT,NPTS,M,NSP)
C00005 00003		SUBROUTINE INVFL(X,NP,M,RO,ERRN)
C00007 00004		SUBROUTINE FFTP(X,Y,M,L)
C00008 00005		SUBROUTINE RBITS(X,Y,M)
C00010 00006		SUBROUTINE USCRM(X,Y,M)
C00011 ENDMK
C⊗;
	SUBROUTINE LPC(AIFFY,SPT,NPTS,M,NSP)
C	THE PARAMETERS IN THE LIST ARE DEFINED AS FOLLOWS:
C	AIFFY←THE INPUT DATA ARRAY OF REAL NUMBERS...A TIME SERIES
C	NPTS←THE NUMBER OF POINTS IN THE TIME SERIES
C	M←   THE NUMBER OF FILTER COEFFICIENTS DESIRED.
C	SPT← THE ARRAY IN WHICH THE REAL SPECTRUM IS RETURNED.
C	NSP← THE NUMBER OF POINTS IN THE REAL SPECTRUM.
	DIMENSION AIFFY(1),CF(50),SPT(1),X(512),Y(512)
	TPI=2.*355./113.
10	NP=NPTS-1
C	DON'T CALCULATE THE DIFFERENCE SIGNAL
	DO 11 J=1,NP
11	X(J)=AIFFY(J)
C 11	X(J)=AIFFY(J+1)-AIFFY(J)
C	APPLY A HAMMING WINDOW TO THE TIME SERIES.
	DT=NP
	X(NPTS)=X(NPTS-1)*.08
		DO 12 J=1,NP
		T=J-1
12		X(J)=X(J)*(.54-.46*COS(TPI*T/DT))
C	CALCULATE AUTOCORRELATION COEFFICIENTS AND SOLVE
C	FOR INVERSE FILTER COEFFICIENTS
100	DO 105 J=1,NPTS
105	Y(J)=X(J)
	CALL INVFL(Y,NPTS,M,RO,ERRN)
C	Y ARRAY RETURNS 1.,A1,A2,...,AM,0.,....0.
	MM=M+1
	DO 125 J=1,MM
125	CF(J)=Y(J)
300	DO 303 J=1,512
	X(J)=0.
303	Y(J)=0.
	MM=M+1
	DO 310 J=1,MM
310	X(J)=CF(J)
C	SHUFFLE DATA FOR CALC OF REAL TRANSFORM
	MM2=MM/2 +1
		DO 320 K=1,MM
		K1=2*K-1
		K2=2*K
		X(K)=X(K1)
320		Y(K)=X(K2)
	NMOD=1
323	IF (2**NMOD.GE.NSP) GO TO 326
	NMOD=NMOD+1
	GO TO 323
326	IF ((2*MM2).LT.MM) MM2=MM2+1
	MOD=1
330	IF (2**MOD.GE.MM2) GO TO 340
	MOD=MOD+1
	GO TO 330
340	CALL FFTP(X,Y,NMOD,MOD)
	CALL RBITS(X,Y,NMOD)
	CALL USCRM(X,Y,NMOD)
C	CALCULATE LOG OF RECIPROCAL INVERSE FILTER SPECTRUM
		DO 350 J=1,NSP
350		X(J)=X(J)*X(J)+Y(J)*Y(J)
		BIG=-1.0E+20
		DO 360 J=1,NSP
		IF (X(J).EQ.0) X(J)=1.*10E-36
		SPT(J)=-10.*ALOG10(X(J))
360		IF(BIG.LT.SPT(J))BIG=SPT(J)
		DO 377 J=1,NSP
377		SPT(J)=SPT(J)-BIG
	RETURN
	END
	SUBROUTINE INVFL(X,NP,M,RO,ERRN)
	DIMENSION X(512)
	DIMENSION F(50),TF(50),A(50),R(50)
C	CALCULATION OF M+1 LENGTH AUTOCORRELATION SEQUENCE
C	FROM THE INVERSE FILTER INPUT SEQUENCE
	MP1=M+1
		DO 11 JJ=1,MP1
		J=JJ-1
		MMJ=NP-J
		SS=0.
			DO 10 I=1,MMJ
			IPJ=I+J
10			SS=SS+X(I)*X(IPJ)
11		R(JJ)=SS
	RO=R(1)
C	FORTRAN IMPLEMENTATION OF LEVINSON'S METHOD FOR
C	SOLVING THE AUTOCORRELATION MATRIX
	F(1)=1.
	ALPHA=R(1)
	BETA=R(2)
	A(1)=-R(2)/R(1)
	GAMMA=A(1)*R(2)
		DO 1 N=2,M
		NM1=N-1
		C=-BETA/ALPHA
		IF (N-2) 2,2,3
3			DO 4 J=2,NM1
			NN=N-J+1
4			TF(J)=F(J)+C*F(NN)
			DO 5 J=2,NM1
5			F(J)=TF(J)
2		F(N)=C*F(1)
		ALPHA=ALPHA+C*BETA
		BETA=0.
			DO 6 J=1,N
			NN=N-J+2
6			BETA=BETA+F(J)*R(NN)
		Q=-(R(N+1)+GAMMA)/ALPHA
			DO 7 J=1,NM1
			NN=N-J+1
7			A(J)=A(J)+Q*F(NN)
		A(N)=Q*F(1)
		GAMMA=0.
			DO 8 J=1,N
			NN=N-J+2
8			GAMMA=GAMMA+A(J)*R(NN)
1		CONTINUE
C	CALCULATE NORMALIZED ERROR ERRN
C	SM=0.
C		DO 77 J=1,M
C 77		SM=SM+A(J)*R(J+1)
C	ERRN=1.+SM/R(1)
C	PLACE COEFFICIENTS IN PROPER LOCATIONS OF THE
C	REAL VECTOR X FOR FFT-ING TO GET INVERSE SPECTRUM
		DO 51 J=1,M
51		X(J+1)=A(J)
	X(1)=1.
	MP1=MP1+1
		DO 123 J=MP1,512
123		X(J)=0.
	RETURN
	END
	SUBROUTINE FFTP(X,Y,M,L)
	DIMENSION X(512),Y(512)
	N=2**M
	L2=2**L
	        DO 1 LO=1,M
	LMX=2**(M-LO)
	LMM=LMX
	LIX=2*LMX
	SCL=6.283185/LIX
C     TEST FOR PRUNING
	IF (LO-M+L) 20,30,30
20	LMM=L2
30	        DO 1 LM=1,LMM
	ARG=(LM-1)*SCL
	C=COS(ARG)
	S=SIN(ARG)
	        DO 1 LI=LIX,N,LIX
	J1=LI-LIX+LM
	J2=J1+LMX
	T1=X(J1)-X(J2)
	T2=Y(J1)-Y(J2)
	X(J1)=X(J1)+X(J2)
	Y(J1)=Y(J1)+Y(J2)
	X(J2)=C*T1+S*T2
1	Y(J2)=C*T2-S*T1
	RETURN
	END
	SUBROUTINE RBITS(X,Y,M)
	DIMENSION X(512),Y(512),L(9)
	EQUIVALENCE(L9,L(1)),(L8,L(2)),(L7,L(3)),(L6,L(4)),
     1 (L5,L(5)),(L4,L(6)),(L3,L(7)),(L2,L(8)),(L1,L(9))
C       PERFORM BIT REVERSAL
	    DO 70 J=1,9
	    L(J)=1
	    IF (J-M) 71,71,70
71	    L(J)=2**(M+1-J)
70	    CONTINUE
	JN=1
	    DO 60 J1=1,L1
	    DO 60 J2=J1,L2,L1
	    DO 60 J3=J2,L3,L2
	    DO 60 J4=J3,L4,L3
	    DO 60 J5=J4,L5,L4
	    DO 60 J6=J5,L6,L5
	    DO 60 J7=J6,L7,L6
	    DO 60 J8=J7,L8,L7
	    DO 60 JR=J8,L9,L8
	IF (JN-JR) 61,61,62
61	R=X(JN)
	X(JN)=X(JR)
	X(JR)=R
	FI=Y(JN)
	Y(JN)=Y(JR)
	Y(JR)=FI
62	JN=JN+1
60	CONTINUE
	RETURN
	END
	SUBROUTINE USCRM(X,Y,M)
	DIMENSION X(512),Y(512)
	I=2**M
	LO2=I/2
	SA=3.141593/I
	X(1)=2.*(X(1)+Y(1))
	Y(1)=0.
	LLL=LO2+1
	X(LLL)=2.*X(LLL)
	Y(LLL)=-2.*Y(LLL)
		DO 33 K=2,LO2
		J=I-K+2
		ANG=SA*(K-1)
		C=COS(ANG)
		S=SIN(ANG)
		AA=X(K)+X(J)
		BB=Y(K)-Y(J)
		CC=X(K)-X(J)
		DD=Y(K)+Y(J)
		XR=C*DD-S*CC
		XI=S*DD+C*CC
		X(K)=AA+XR
		Y(K)=BB-XI
		X(J)=AA-XR
		Y(J)=-BB-XI
33	CONTINUE
	RETURN
	END